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Abstract 

Most  variational  models  for  multi-phase  image  segmentation  are 
non-convex  and  possess  multiple  local  minima,  which  makes  solving 
for  a  global  solution  an  extremely  difficult  task.  In  this  work,  we  pro¬ 
vide  a  method  for  computing  a  global  solution  for  the  (non-convex) 
multi-phase  piecewise  constant  Mumford-Shah  (spatially  continuous 
Potts)  image  segmentation  problem.  Our  approach  is  based  on  using 
a  specific  representation  of  the  problem  due  to  Lie  et  al.  [27].  We 
then  rewrite  this  representation  using  the  dual  formulation  for  total 
variation  so  that  a  variational  convexification  technique  due  to  Pock 
et  al.  [30]  may  be  employed.  Unlike  some  recent  methods  in  this  di¬ 
rection,  our  method  can  guarantee  that  a  global  solution  is  obtained. 

We  believe  our  method  to  be  the  first  in  the  literature  that  can  make 
this  claim.  Once  we  have  the  convex  optimization  problem,  we  give 
an  algorithm  to  compute  a  global  solution.  We  demonstrate  our  algo¬ 
rithm  on  several  multi-phase  image  segmentation  examples,  including 
a  medical  imaging  application. 

*  Authors  are  with  the  Department  of  Mathematics,  University  of  California,  Los  Ange¬ 
les.  Email:  {ethan, chan,xbresson}@math.ucla. edu.  This  work  was  supported  in  part 
by  ONR  N00014-03-1-0071  and  ONR  MURI  subcontract  from  Stanford  University. 
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1  Introduction 


Image  segmentation  is  one  of  the  most  fundamental  problems  in  computer 
vision.  The  task  is  to  assign  a  label  to  each  point  in  a  given  image  (i.e., 
partition  the  image  into  phases  or  segments)  such  that  the  labeling  is  optimal 
with  respect  to  a  particular  model.  It  is  difficult  to  compute  a  globally 
optimal  solution  because  the  models  often  lead  to  non-convex  variational 
problems.  Most  techniques  instead  devise  algorithms  to  efficiently  compute 
solutions  which  are  only  locally  optimal. 

The  labels  assigned  during  image  segmentation  can  be  defined  either  on 
a  discrete  grid  or  in  a  spatially  continuous  domain.  In  the  discrete  setting, 
interest  for  guaranteed  global  solutions  of  vision  problems  has  grown  ever 
since  the  stochastic  approximate  methods  in  the  seminal  work  of  Geman  and 
Geman  [21].  Greig  et  al.  [22]  were  the  first  to  discover  that  the  classical 
Ford-Fulkerson  min- cut /max- flow  correspondence  [20]  could  provide  global 
algorithms  when  the  set  of  labels  is  binary  and  the  model  satisfies  certain 
criteria.  The  work  of  [23]  showed  that  a  class  of  multi-label  problems  can 
be  globally  solved,  but  this  does  not  include  one  of  the  most  celebrated 
image  segmentation  models:  the  piecewise  constant  Mumford-Shah  model 
[28]  (known  as  the  Potts  model  [31]  for  the  discrete  problem).  Recently,  for 
this  model,  Bae  and  Tai  [2]  provide  a  global  graph  cut  method  for  the  case 
of  four  labels  under  some  technical  assumptions.  In  this  paper,  however,  our 
focus  is  the  continuous  setting;  see  [24]  for  some  discussion  on  the  comparison 
between  discrete  and  continuous  methods  in  optimization  problems. 

For  the  two-phase  case  (i.e.,  assigning  one  of  just  two  labels,  foreground 
or  background,  to  each  point  in  the  continuous  domain)  Chan  et  al.  [12] 
developed  an  equivalent  convex  formulation,  thereby  constructing  an  algo¬ 
rithm  that  guarantees  to  find  a  globally  optimal  solution.  The  method  relics 
on  a  relaxation  of  the  binary  labeling  function  which  remarkably  still  yields 
a  solution  that  is  binary.  However,  this  technique  is  not  directly  applicable 
to  the  multi-phase  case,  which  has  remained  an  open  problem. 

Inspired  by  [12],  there  has  been  significant  research  devoted  to  globally 
solving  the  multi-phase  segmentation  problem  (or  equivalent  multi-labeling 
problems).  Zach  et  al.  [35]  proposed  a  relaxation  approach  along  with  a 
decoupling  to  yield  a  convex  optimization  problem,  but  could  not  guarantee 
a  global  optimum  of  the  original  problem.  Lellmann  et  al.  [26]  used  a  similar 
method,  with  a  slightly  different  regularization  and  splitting  technique,  but 
this  had  comparable  limitations.  Most  recently,  Pock  et  al.  [29]  used  the  dual 
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formulation  of  total  variation  (TV),  enforcing  a  convex  constraint  on  the  dual 
variables.  Their  relaxation  strictly  dominates  these  previous  approaches,  in 
the  sense  that  the  set  of  functions  over  which  their  convex  minimization  is 
conducted  contains  that  of  the  earlier  methods,  which  results  in  a  better 
approximation.  Nevertheless,  this  method  can  still  not  guarantee  a  globally 
optimal  solution. 

The  goal  of  this  paper  is  to  provide  a  convex  formulation  and  an  algorithm 
that  is  guaranteed  to  compute  a  global  solution  for  the  multi-phase  piecewise 
constant  Mumford-Shah  model.  To  the  best  of  our  knowledge,  this  paper 
provides  the  first  approach  that  can  make  this  guarantee. 

Our  method  relies  primarily  on  two  existing  approaches  for  image  segmen¬ 
tation  and  computer  vision  problems.  First,  we  use  the  so-called  piecewise 
constant  level  set  method  (PCLSM)  framework  of  Lie  et  al.  [27]  to  formu¬ 
late  the  problem,  which  provides  a  convenient  representation  to  enforce  the 
constraints  in  the  optimization.  Second,  we  use  the  novel  convexification  ap¬ 
proach  of  Pock  et  al.  [30]  (which  was  motivated  by  the  work  of  Ishikawa  in  the 
discrete  Markov  Random  Field  setting  [23])  that  converts  certain  non-co nvex 
variational  problems  to  equivalent  convex  problems  in  a  higher  dimensional 
space.  LIsing  the  dual  formulation  of  TV  written  in  the  PCLSM  framework 
allows  our  method  to  obtain  a  convex-concave  min-max  problem  from  which 
a  global  solution  may  be  computed. 

The  paper  is  organized  as  follows.  In  Section  2,  we  formally  introduce 
the  segmentation  problem  and  explain  the  solution  of  [12]  for  the  two-phase 
case.  In  Section  3,  we  begin  by  describing  the  representation  we  use  for  the 
multi-phase  problem  as  well  as  the  convexification  technique  of  [30] .  Later  in 
the  section,  we  show  our  main  result:  that  we  can  find  a  global  solution  for 
multi-phase  segmentation;  we  also  describe  an  algorithm  for  computing  such 
a  solution.  Section  4  explains  the  numerical  implementation  of  the  algorithm 
and  shows  some  experimental  results.  Finally,  in  Section  5,  we  compare  our 
method  to  previous  work  in  this  direction  and  give  some  concluding  remarks. 
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2  Segmentation  problem  and  two-phase  solu¬ 
tion 


Let  I :  Q  — >  R  be  a  given  image  with  bounded  image  domain  C 
Mumford  and  Shah  [28]  proposed  the  minimization  problem 


inf 

f,c 


(I(x)  —  f(x))2  dx  +  v 


>n\c 


| V/(x)|2  dx  +  p\C\ 


to  hnd  an  optimal  piecewise  smooth  approximation  of  a  general  function 
(the  image  I),  where  v  and  /i  are  fixed  parameters.  The  edge  set  C  C  hi  is  a 
closed  set  that  defines  a  partition  hi  =  such  that  the  restrictions  ft  of 
the  function  /  to  the  segments  flj  are  differentiable.  In  the  simplest  form  of 
the  model,  the  function  /  is  taken  to  be  constant  on  each  segment,  reducing 
to  the  piecewise  constant  Mumford-Shah  problem 

inf  i  f  (c*  —  I(x))2  dx  +  fi \C\ 

Ci'c  {  i  •/<>. 


where  /  =  JT  Cjln;.  This  is  also  known  as  the  Potts  model  [31]  in  the  discrete 
setting,  originating  in  solid  state  physics,  and  hence  we  will  use  the  names 
of  these  models  interchangeably.  The  important  property  of  the  model  is 
imposing  a  regularity  measure  that  favors  the  labels  of  neighboring  points  to 
be  identical  but  interprets  the  value  of  each  label  to  be  immaterial. 

We  will  herein  assume  that  the  optimal  constants  ct  are  known  a  priori 
and  the  number  of  segments  (say,  k  +  1)  is  fixed.  Indeed,  even  making  this 
stringent  assumption  leaves  us  with  the  difficult  non-convex  problem 

k  i  r 

V'-|<9rb|+  /  (Ci  -  I(x))2dx 
i= o  1  J 

where  it  is  implicit  that  f!  =  Ujflj  and  the  flj  are  pairwise  disjoint.  In  this 
work,  we  will  not  consider  the  even  more  challenging  problem  where  either  the 
constants  or  number  of  phases  are  unknown.  Likewise,  the  related  methods 
for  the  multi-phase  problem  mentioned  in  the  Introduction  and  discussed 
further  in  Section  5  also  make  these  reasonable  assumptions. 

To  slightly  generalize  our  method  to  data  terms  other  than  the  one  in  the 
Mumford-Shah  functional,  the  multi-phase  problem  we  will  solve  is 

+  /* 

J  Q,i 


inf 

n0,...,nfc 


\dVti 
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For  example,  when  the  region  descriptors  gi(x)  are  equal  to  (ct  —  I(x ))2 
(where  we  have  absorbed  constant  factors),  we  have  the  case  of  Mumford- 
Shah.  This  generalization  does  not  add  any  complexity  to  our  problem.  Our 
assumption  concerning  the  c*  translates  to  the  assumption  that  the  gi(x)  are 
known  beforehand. 

When  k  —  1,  expression  (1)  above  becomes  the  two-phase  problem 

inf  ( |<9O0|  +  [  g0(x)dx  +  [  g^x)  dxX  .  (2) 

n»cn  l  Jci0  Jn\n0  J 

Many  popular  methods  to  solve  (2)  are  based  on  the  level  set  method  (e.g., 
the  Chan-Vese  method  [14]),  but  there  is  no  guarantee  to  find  a  global  solu¬ 
tion.  In  [12],  Chan  et  al.  proposed  a  method  for  finding  global  minimizers 
through  a  convex  formulation  of  the  problem.  They  showed  that  (2)  can  be 
written  equivalently 

min  {  E(9)  :=  (  |V0|  dx  +  [  (1  —  9(x))g0(x)  +  9(x)gi(x)  dx\ ,  (3) 

^e{o,i}  (  Jn  Jn  J 

where  6 :  — >  {0, 1}  is  a  binary  function  that  defines  the  segmentation: 
x  G  O0  if  9(x)  =  0  and  x  G  fR  if  6{x)  =  1.  However,  even  though  E{6)  is 
convex,  the  minimization  is  done  over  a  non-convex  set  of  binary  functions. 

The  important  final  step  is  that  a  relaxation  may  be  taken  to  allow  9  € 
[0, 1]  without  changing  the  minimum.  Indeed,  letting  6*  be  any  minimizer  of 
the  binary  problem  (3),  then  any  9  G  [0, 1]  satisfies  E{9*)  <  E{9).  Moreover, 
a  solution  of  (3)  can  be  obtained  by  finding  a  solution  of  the  relaxed  problem 
and  then  thresholding.  Namely,  letting  now  9*  be  a  minimizer  of  the  relaxed 
problem,  the  binary  function  l{e»*>t}  is  a  minimizer  of  (3)  for  any  t  G  (0,1). 
This  relaxation  completes  the  conversion  of  the  original  problem  (2)  to  a 
convex  minimization  problem.  We  say  the  relaxation  is  exact  because  a  true 
solution  of  the  original  problem  can  be  exactly  recovered  from  the  relaxed 
problem. 

Relaxation  in  this  context  dates  back  to  Strang’s  work  [33]  on  maximal 
flows.  The  critical  property  for  the  functional  of  a  problem 

min  Fix) 

a;G{0,l} 

to  possess  in  order  to  have  an  exact  relaxation  to  x  G  [0, 1]  is  a  “generalized 
co-area  formula”  of  the  form 

F(x)  =  [  F(l{x>t})  dt. 

Jo 
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The  article  [7]  contains  a  nice  presentation  of  the  proof  that  the  relaxation  is 
exact.  As  it  turns  out,  the  fact  that  the  relaxation  does  not  have  a  straight¬ 
forward  extension  to  vector-valued  functions  is  the  obstacle  in  many  of  the 
recent  multi-phase  approaches  mentioned  in  the  introduction;  see  Section  5 
for  further  discussion. 


3  Extension  to  multi-phase 


3.1  Multi-phase  problem  representation 

We  return  now  to  developing  the  tools  to  globally  solve  the  multi-phase 
problem 

+  [  9i 

Moving  forward,  it  becomes  advantageous  to  represent  this  problem  as  an 
equivalent  optimization  problem  over  functions  rather  than  over  partitions 
of  12.  In  the  literature,  there  are  a  number  of  representations  of  the  segmen¬ 
tation  problem  as  such,  e.g.  [34,  15].  For  reasons  which  we  hope  to  clarify 
throughout  this  paper,  we  use  the  PCLSM  representation  of  Lie  et  al  [27]  in 
order  to  globally  solve  the  problem.  We  represent  the  partition  12  =  U*12j  by 
a  piecewise  constant  labeling  function  u :  12  — >  (0, . . . ,  k}  with  the  property 
that  u  =  i  on  12*.  For  convenience,  we  introduce  the  notation  ipi(u)  :=  l{u=i} 
for  the  k  +  1  characteristic  functions  induced  by  u.  Clearly  the  function  u  is 
given  by  u  =  Ei=i  b/y •  The  characteristic  functions  can  further  be  utilized 
to  express  the  boundary  lengths  of  the  segments: 


inf 

^0  )••* 


]T|<9Lh 


i= 0 


1(912*1  =  /  |V^(w)l  dx. 

Jn 

This  is  a  consequence  of  the  fact  that  the  TV  of  the  characteristic  function  of 
a  set  is  its  perimeter.  Under  this  representation,  the  segmentation  problem 
(4)  becomes  the  constrained  optimization  problem 


mm 
£7— >{0 


u)  I  +  tjji(u)gi(x)  dx 


Next,  because  we  eventually  wish  to  obtain  a  continuous  convex  formulation, 
we  relax  and  let  u  G  M.  To  ensure  that  the  function  u  takes  values  in 
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{0, . . . ,  k},  we  introduce 


k 

K(u)  =  JJ(u-i), 

i= 0 

so  that  if  a  function  u:  Q  — >  M  satisfies  identically  K{u)  =  0,  then  there 
exists  a  unique  i  G  {0 . . . ,  k}  for  every  x  G  O  such  that  u(x)  =  i.  We  also 
use  the  interpolation  formulas 


(u~j) 
(i  -  j)  ’ 


which  coincide  with  the  characteristic  functions  on  the  constraint  set.  Our 
problem  is  now 

min  \y  [  \V^i(u)\+ ^i(u)gi(x)dx\  s.t.  K(u)  —  0.  (5) 

u:  o^[oy]  Jn  J 

Notice  there  is  no  harm  replacing  the  condition  mGI  with  u  G  [0,  k]  because 
the  constraint  K(u)  already  demands  u  G  {0, . . . ,  k}.  Globally  solving  this 
problem  appears  to  be  quite  formidable  since  both  the  objective  function 
and  the  feasible  set  are  non-convex.  The  remedy  is  a  technique  for  convert¬ 
ing  certain  non-convex  problems  to  convex  ones,  which  we  describe  in  the 
following  subsection. 


3.2  Convex  relaxation  by  functional  lifting 

We  briefly  review  the  method  of  [30]  to  solve  minimization  problems  of  the 
form 


min  /  |  Vu|  +  p(x,  u(x))  dx,  (6) 

u-.  n->r:=[7min,7max]  Jn 

where  p:  x  T  — >  M  may  be  non-convex.  The  method  transforms  this 
non-convex  problem  into  an  equivalent  convex  problem  through  a  change 
of  independent  variables  and  a  relaxation  on  the  constraints.  While  our 
problem  (5)  is  not  quite  of  this  form,  we  will  eventually  use  this  method 
in  our  technique.  In  our  application,  we  will  frequently  take  7min  =  0  and 
7max  =  k  since  these  are  the  endpoints  of  the  feasible  set  of  (5). 

Simply  put,  the  idea  is  to  reformulate  the  problem  in  terms  of  the  super¬ 
level  set  function  0:  x  T  — >  {0, 1}  defined  by  0(x,7)  =  1{u(j,)>7}(x)  in  such 
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a  way  that  the  objective  function  of  the  reformulated  problem  is  convex.  The 
function  u  may  then  be  recovered  from  0  from  the  layer-cake  formula 


u(x)  =  7min+  /  00c,  7)^7 . 


(7) 


One  tool  used  to  accomplish  this  reformulation  is  the  co-area  formula  [18,  19] 


\Wu\dx 


|  Vcpld'ydx. 


The  other  key  observation  is  that  S(u(x)  —  7)  =  |<970(:r,  7)].  Indeed,  since  4 > 
is  the  super-level  set  function  of  u, 


<t>(l,x) 


1  if  7  <  u(x), 
0  otherwise. 


Hence,  the  derivative  of  0  with  respect  to  7  is  zero  except  when  u(x)  =  7, 
where  a  jump  occurs  and  the  magnitude  of  this  derivative  is  a  delta  function 
centered  at  u(x)  —  7.  It  follows  that 


p(x,u(x))dx  —  /  /  p(x,  7)  |<970|  d'y  dx. 


q  Jr 


Thus,  when  0  is  the  super-level  set  function  of  u,  we  have 


/  \Vu\  +  p(x,  u(x))  dx  =  /  /  |V0|  +  p(x,  7)  |<970|  d'ydx.  (8) 

Jn  J n  Jr 

This  implies  that  we  may  solve  (6)  from 

min  /  /  |V0|  +  p(x,7)  |970|  d'ydx,  (9) 

<^{0,1}  Jn  J r 

as  long  as  we  require  0  to  be  a  super-level  set  function,  i.e.,  0(x,  7m;n)  = 
l,0(a:,7max)  =  0,  and  <p(x,  •)  is  non- increasing.  We  then  use  (7)  to  obtain  u 
from  0.  Observe  that  the  objective  function  of  (9)  is  convex  in  0  because  the 
non-convex  function  p  does  not  depend  on  0;  the  non-convexity  has  somehow 
been  integrated  out. 

The  authors  of  [30]  call  this  technique  functional  lifting  because  we  have 
increased  the  dimension  of  the  arguments  of  the  functions  over  which  the 


minimization  is  taken.  Indeed,  rather  than  minimize  over  functions  u(x), 
we  instead  minimize  over  functions  0(:r,  7).  The  idea  was  motivated  by  the 
discrete  approach  of  Ishikawa  [23]  in  which  an  appropriate  auxiliary  graph 
with  an  extended  node  set  was  created,  and  standard  (binary)  graph-cut 
computations  (see,  e.g.,  [25])  were  used  to  solve  the  problem. 

Finally,  analogous  to  the  procedure  in  [12],  a  relaxation  on  (p  is  taken  to 
allow  <p  £  [0, 1].  This  is  crucial  to  ultimately  obtain  a  convex  minimization 
problem.  The  form  of  the  integrand  of  (9)  allows  for  essentially  the  same 
proof  as  that  of  [12]  to  show  the  relaxation  is  exact,  in  the  sense  described  at 
the  end  of  Section  2.  In  other  words,  thresholding  a  minimize!'  of  the  convex 
problem 

min  /  /  |V0|  +p(x,7)  |d70|  d'ydx, 

Jn  Jr 

where 

D  =  {(f>:  Q.  x  T  ->  [0, 1]:  <p(x,  7min)  =  1,  (p(x,  7max)  =  0}, 

yields  a  minimize!'  of  (9),  which  in  turn  via  (7)  results  in  a  solution  of  the 
original  non-convex  problem.  It  is  worth  mentioning  that  because  the  result¬ 
ing  objective  function  is  convex  but  not  strictly  so,  there  is  in  general  not  a 
unique  minimizer. 

Despite  the  effectiveness  of  this  method,  it  does  not  directly  apply  to 
our  segmentation  model.  The  first  term  in  (6)  is  the  TV  of  the  labeling 
function  u.  This  penalizes  boundaries  proportional  to  the  difference  between 
the  values  in  the  neighboring  segments.  Instead,  our  model  (as  well  as  many 
others)  wants  to  penalize  all  boundaries  equally  since  the  value  of  the  labels 
should  be  arbitrary. 

3.3  Convex  method  for  multi-phase  problem 

We  are  now  ready  to  explain  our  method  for  finding  a  global  solution  of 
(5).  To  enforce  the  constraint  K{u)  =  0,  we  use  an  augmented  Lagrangian 
method  [3].  Using  results  from  the  theory  of  constrained  optimization,  if  we 
find  a  sequence  {uj}  of  global  minimizers  of 


min 

u£[0,fc] 


(10) 
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for  appropriate  choices  of  the  multipliers  A j  €  L2(fl)  and  the  penalty  param¬ 
eters  Tj  >  0,  then  any  limit  point  of  the  sequence  will  be  a  global  minimizer 
of  (5).  For  example,  one  condition  that  guarantees  such  a  global  minimizer 
is  when  the  set  of  multipliers  is  bounded  and  the  penalty  parameters  satisfy 
ry  <  r2  <  ■  ■  ■  and  Tj  — >•  oo  [3,  Prop.  2.1].  In  practice,  appropriately  chosen 
values  for  Xj ,  ry  can  provide  better  convergence  results;  see  Section  4.  It 
remains  to  find  a  global  minimizer  of  (10)  for  fixed  j. 

To  this  end,  we  use  the  dual  formulation  of  TV  [13,  8,  9].  We  introduce 
dual  variables  pi :  Q  — >  M.d  (recall  that  U  C  Md)  for  i  —  0, . . . ,  k  so  that  (10) 
becomes 


min  max 
«e[0,fc]  |pi]<l 


where  we  have  put 


y]  (—  di v pi  ■  ipi(u))  +  Wj{x ,  u{x))  dx 


i= o 


(ii) 


k 

Wj{x,  7)  =  E  +  \K( 7)  +  ^A'(y)2. 

i=0 

The  norm  in  the  maximization  constraint  in  the  above  expression  refers  to 
the  L2  norm  on  the  dual  variables. 

Using  the  functional  lifting  method  described  in  Section  3.2,  we  see  that 
if  0  is  the  super-level  set  function  of  u,  then 


„  k 

max  /  y  (—  divpj  •  ipi(u))  +  Wj(x,  u{x ))  dx 


is  equal  to 


max 
bil<i  I  Jn  Jo 


E  (“  divP* '  ^(7))  +  Wj(x,  7)  |<970|  dy  dx  l , 


vi=0 


and  hence  (11)  is  equivalent  to 


min  max  < 

<j>&D  |Pi|<l 


k  /  k 


a  Jo 


E  (~  divP»  •  V,i( 7))  +  wj(x ,  7)  1^701  dldx 


.i=0 


UO.p) 
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where 

D  =  {(f>:  fl  x  T  — >  [0, 1] :  cj)(x,  0)  =  1,  <f{x ,  k )  =  0}.  (12) 


Put 


X  =  {p  =  (p0,...,pk):  n^Mdx(fc+1):  \Pi\  <  1  Vi  =  0,  ...,&}  .  (13) 


We  have  an  optimization  problem  over  the  set  D  x  X,  where  D  and  X  are 
compact,  convex  subsets  of  linear  topological  spaces,  such  that  $j(0,  •)  is 
concave  on  X  for  all  0  G  D  and  <!>.,'(•,  p)  is  convex  on  D  for  all  p  e  X.  By 
Sion’s  minimax  theorem  [32],  the  minimization  and  maximization  operations 
may  be  freely  interchanged. 

We  summarize  our  main  result  in  the  following  theorem,  whose  proof 
follows  from  the  entirety  of  this  section. 


Theorem  1.  The  multi-phase  segmentation  problem 


equivalently 


mm 


u :  fi— d0,fc]  |  —  J  fi 


£/  |V0i(w)|  +  4>i(u)gi(x)  dx  >  s.t.  K(u)  =  0, 


can  be  globally  solved  by  finding  a  limit  point  u*  of  a  sequence  of  functions  u.j, 
constructed  as  follows.  Suppose  ( { A ^ } ,  {r:)})°V1  is  such  that  {Aj}  is  bounded, 
r'i  <  r2  <  ■  ■  ■ ,  and  r,  — >•  oo.  For  each  j ,  globally  solve  the  convex  optimization 
problem 

min  max  $^(0,  p) 

<fi£D  pSX 

to  obtain  a  function  0*  and  then  form  the  binary  function  f>j  :=  1  {</,*>*}  for 
any  t  G  (0, 1).  Finally,  compute  the  solution  Uj  via 

nk 

Uj(x)  =  /  0j(rr,7 )drr- 
Jo 
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3.4  Algorithm 


Let  us  now  describe  an  algorithm  to  obtain  a  global  optimal  solution  of  (1), 
as  guaranteed  by  Theorem  1.  Our  objective  is  simply  to  describe  a  global 
algorithm  resulting  from  our  theory  rather  than  developing  the  most  efficient 
algorithm  possible,  which  we  leave  for  future  work. 

For  the  inner  maximization  step,  observe  that  the  terms  in  our  primal- 
dual  objective  function  dependent  on  the  dual  variables  pi  decouple.  In  other 
words,  to  conduct  the  maximization,  we  may  separately  solve  for  each  v. 


max 

N<i 


div pi  ■  ipi( 7)  |<970|  djdx>  , 


or  equivalently, 


max 

|p*l<i 


Ail)  dj 


This  may  be  solved  using  the  iterative  scheme  [10]: 


P 


n+1 


=  Vx[pni+TpV 


Ail)  \d,r\ d 7 


(14) 


where  rv  denotes  the  step  size  of  updates  for  the  pi  variables.  The  operator 
Vx  is  the  projection  onto  the  set  X  defined  in  (13),  i.e. , 


Vx(ci) 


qo 


max(q0, 1)  ’  "  '  ’  max(qfc,  1) 


for  all  q:  — >  Mdx(fc+1b 

For  the  outer  minimization  problem,  we  use  an  explicit  gradient  descent 
method  derived  from  the  corresponding  Euler-Lagrange  equation.  Again, 
let  us  emphasize  that  this  is  just  one  method  to  perform  this  minimization 
subproblem  and  certainly  not  the  most  efficient.  It  was  instead  chosen  for 
ease  of  implementation.  Put 


k 

Wj(x,  7,p)  =  ^(-divpj(x)  •ipi('y))+wj(x,'y). 

i= 0 


The  problem  is  then,  for  fixed  p, 


mm 

<f>£D 


W(x,  7,  p)  |<97</>!  dr)  dx 


n  Jo 


) 
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Introducing  an  artificial  time  t,  the  gradient  descent  PDE  is 


dp 

dt 


-V7- 


W{x,  7, 


We  write  the  divergence  operator  V7-  =  divT  instead  of  d1  since  these  will  be 
different  in  the  finite  difference  numerical  scheme.  To  avoid  the  degenerate 
case  when  d^p  =  0,  we  regularize  this  equation  by  some  small  e  >  0  to  yield 


dp 

dt 


-V, 


W(x,  7,p) 


<970 


e2  +  (<970)2 


This  gives  the  iterative  scheme 


r+1  =  Pg{r+T,vT.(ty(l,7,p)_||_)},  (15) 

where  Vd  is  the  Euclidean  projection  onto  the  convex  set  D  defined  in  (12). 
Summarizing,  we  propose  the  following  algorithm. 


Algorithm  1.  1.  Fix  step  sizes  rp >  0,  convergence  tolerances  0  < 

eu,e^  <  1,  and  regularization  parameter  0<e€l. 


2.  Let  j  be  the  iteration  number  for  the  functio7i  u,  multiplier  A,  and 
penalty  parameter  r.  Initialize  u°,  A0,  and  r°.  Until  convergence 
\'u]+]  —  u3 :|  <  eu,  update  u3 ,  X3 ,  and  r3  with  the  following: 


(a)  Let  n  be  the  iteration  number  for  the  variables  p  and  pt  (i  = 
0, . . .  ,k).  Initialize  p°  as  the  super-level  set  function  of  u3 .  Ini¬ 
tialize  =  0.  Until  convergence  \ pn+1  —  pn\  <  update  pn  and 
pf  with  the  following: 


(b) 


i.  p^step:  For  each  i,  update  using  (14). 

ii.  p-step:  Update  using  (15). 

Threshold  pn+1  to  then  obtain  u3+l  from  the  layer-cake 

formula 

rk 


uJ+1  =  I  dj. 


Choose  AJ+1  and  r3+x  to  satisfy  the  conditions  of  the  theorem. 


Note  that  the  Euclidean  projection  Vd  can  be  implemented  by  a  simple 
point-wise  truncation  operation. 
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4  Numerical  implementation  and  experimen¬ 
tal  results 

To  simplify  the  discussion  of  the  numerical  implementation  of  the  algorithm, 
we  consider  only  the  d  —  2  case  and  let  =  (0,  l)2.  We  use  the  standard 
discretization  in  the  spatial  domain 

Qh  =  {l,...,N}x{l,...,N} 

with  spatial  step  size  h  =  1/N,  and  use  the  discretization 

rA^  =  {(^-l)A7:  l<g<N,} 

with  iV7  =  1  +  k/ Ay.  A  spatially  continuous  function  /:  0  x  T  — >  M  is  now 
approximated  by  a  discrete  function  fh'Ary,  from  which  we  will  often  omit  the 
superscripts.  For  /  :  {1, . . . ,  N}2  x  {1, . . . ,  1V7}  — >  M,  we  use  the  notation 
:=  fihJid)-  We  have  the  following  finite  difference  formulas  for  the 
derivative  operators: 


(div7  f)i,j,g  ( fi,j,g 
(plf)i,j,g  =  (fi,j,g+ 1  — 

tV/)!,,7  =  ( ./V  - 1  ../.</  - 

,(V/)fJ>9  =  (fi,j+l,g  -  fi,j,g)/h 


with  appropriate  boundary  conditions  (see  [9]).  In  our  experiments,  we  are 
given  a  discrete  image  defined  on  Qh  and  let  h  —  1.  We  set  Ay  =  0.25. 

We  mentioned  in  the  previous  section  that  the  augmented  Lagrangian 
variables  A-7  and  r3  should  satisfy  the  restrictions  of  Theorem  1  in  order  to 
guarantee  that  the  algorithm  obtains  a  globally  optimal  solution.  In  practice, 
these  variables  may  be  chosen  differently  to  increase  the  efficiency  of  the 
algorithm.  For  example,  in  the  augmented  Lagrangian  technique  of  [27],  the 
penalty  parameters  r3  are  taken  to  be  equal  to  some  constant  r  for  all  j,  and 
the  Lagrange  multipliers  are  updated  as  A-7"1-1  =  \3  +  rK (u3)]  we  elect  to  use 
this  update  in  our  implementation.  We  let  r  —  1  and  initialize  A0  =  1. 

Our  first  two  experiments  are  given  to  illustrate  that  our  method  can 
obtain  a  global  solution  in  situations  where  previous  methods  may  not.  The 
first  is  well-known  in  the  discrete  graph-cut  literature.  In  [4],  an  experimental 
comparison  of  min-flow/max-cut  algorithms  is  given.  In  particular,  for  the 
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(a)  Input 


(b)  Result 


Figure  1:  Segmentation  of  the  diamond  image. 

multi-labeling  problem  of  the  diamond  picture  in  Figure  1(a),  the  authors  of 
[4]  explain  that  the  problem  is  NP-hard  and  instead  approximate  the  solution 
using  the  iterative  expansion  method  of  [5].  The  method  is  proved  to  find 
an  approximate  solution  within  a  factor  of  two  of  the  optimal  value  of  the 
objective  function,  but  cannot  find  a  global  solution.  On  the  other  hand,  our 
method  produces  the  global  minimum.  The  segmentation  result  is  given  in 
Figure  1(b). 

Next,  we  consider  a  synthetic  three-phase  segmentation  problem  whose 
optimal  solution  is  given  by  a  triple  junction.  This  example  is  motivated  by 
an  example  in  [30] .  The  authors  of  [30]  demonstrated  that  their  method  finds 
a  solution  that  is  closer  to  the  optimal  solution  than  the  previous  techniques 
of  [35,  26].  However,  the  solution  obtained  there  is  not  binary  (i.e. ,  it  does  not 
uniquely  assign  exactly  one  label  to  each  point)  and  hence  is  not  a  solution 
of  the  original  problem.  In  contrast,  our  method  produces  a  valid  solution 
by  construction  and  our  theory  guarantees  its  optimality.  The  segmentation 
result  is  shown  in  Figure  2. 

We  now  illustrate  our  algorithm  on  two  real-world  examples.  Figure  3(a) 
shows  an  MRI  brain  image  common  to  medical  imaging  applications.  Fig¬ 
ure  3(b)  shows  the  result  with  four  phases.  We  show  each  phase  in  a  Fig¬ 
ures  3(c)-(f).  Finally,  we  apply  the  algorithm  to  a  color  image  in  Figure  4 
using  four  phases. 


(a)  Input 


(b)  Result 


Figure  2:  Segmentation  of  a  synthetic  image.  The  optimal  solution  is  given 
by  a  120°  triple  junction. 


Recall  that  our  method  assumes  that  the  region  descriptors  (which,  in 
the  case  of  piecewise  constant  Mumford-Shah,  are  characterized  by  the  mean 
intensities  of  each  segment)  are  known  a  priori.  In  our  experiments,  we 
apply  a  fc- means  clustering  algorithm  to  determine  these  values.  Variations 
of  our  algorithm  could  include  updating  these  values  during  the  segmentation 
process,  but  at  the  moment  our  theory  would  not  guarantee  an  optimal 
solution.  Another  alternative  would  be  to  run  the  entire  algorithm  several 
times,  updating  these  values  after  convergence,  to  yield  slightly  better  results. 

5  Comparison  with  other  approaches 

In  the  previous  sections  we  described  a  method  for  globally  solving  the 
multi-phase  segmentation  problem  under  the  Mumford-Shah/Potts  model. 
Our  method  relied  upon  a  particular  way  of  representing  the  optimization, 
namely,  the  PCLSM  representation  of  Lie  et  al.  [27].  We  now  contrast  our 
approach  with  the  other  ones  mentioned  in  the  Introduction,  none  of  which 
can  guarantee  a  global  solution. 
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(a)  Input 


(b)  Result 


(c)  Phase  1 


(d)  Phase  2 


(e)  Phase  3 


(f)  Phase  4 


Figure  3:  Segmentation  of  an  MRI  brain  image  into  four  phases. 
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(a)  Input  (b)  Result 

Figure  4:  Segmentation  of  a  natural  color  image  into  fonr  phases. 


5.1  Convex  relaxation  over  simplex  (Zach  et  al.  [35], 
Lellmann  et  al.  [26]) 

Recall  that  we  decompose  the  labeling  function  u:  12  — ►  {0, . . .  ,k}  with 


A  = 


1  if  u  —  i, 

0  otherwise. 


This  means  that  u  =  JT  iA-  Rather  than  proceed  as  in  Section  3.3,  consider 
instead  the  formulation 


'p=(V'o,— 


min  |f(T):=V  /  \VA\  +  A9i(x)  dx  \  , 
>o,...,^k)es  [  Jn  J 


(16) 


where 


s=  |(/o,...,/fc):  n  —  (0,  l}*-*-1 :  =  1  for  all  x  G  fij  . 

The  natural  multi-dimensional  analogue  of  the  one- dimensional  convex  re¬ 
laxation  approach  (like  that  of  [12]  or  [30])  would  be  to  minimize  over  the 
convex  probability  simplex 

C  =  j  C/b,  ■  ■  ■ ,  fk) :  ft  -  [0,  l]fc+1  :'%2fi(x)  =  1  for  all  xeQ 
l  i=0 
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so  that  we  have  the  convex  minimization  problem 


min 


k  „ 

V  /  |V^i(x)|  +ipi(x)gi(x)dx 
i= o  Jn 


(17) 


This  is  essentially  the  method  proposed  in  [35].  To  solve  (17),  additional 
variables  (v0, . . . ,  Vk) :  hi  — >  Rk  were  introduced  to  decouple  the  regularization 
and  data  terms  [1,  6],  yielding 


min 

C *l>o,-,ipk),(vo,-,vk)£C 


+  X!  Afji  dx 


5 


where  g  >  0  is  a  parameter  that  controls  the  quadratic  approximation  of 
ifii  and  Vi  for  each  i.  Similarly,  Lellmann  et  al.  [26]  considered  the  relaxed 
problem 


Up  to  constant  factors,  this  is  equivalent  to  (17)  when  (^0,  •  •  • ,  V’fc)  £  S.  To 
solve  the  convex  optimization  problem  (18),  a  Douglas- Rachford  splitting 
algorithm  [16]  was  used. 

The  issue  with  both  of  these  methods  is  that  the  minimizers  of  (17) 
and  (18)  cannot  be  guaranteed  to  lie  in  the  set  S.  More  precisely,  if  T*  = 
(^o,  •  •  • ,  ipl)  a  niinimizer  of,  for  example,  (17),  then  it  could  be  the  case 
F(T*)  is  strictly  less  than  the  minimum  of  F  over  S.  It  is  unfortunately  not 
possible  (in  contrast  to  [12]  or  [30])  to  obtain  a  niinimizer  of  (16)  simply  by 
constructing  (l{^»>q, . . . ,  l{^,*>t|)  for  any  t  G  (0,1),  since  this  thresholded 
vector  need  not  even  be  in  S.  It  is  also  not  possible  to  obtain  the  niinimizer 
by  assigning,  for  each  i£fl,  the  unit  vector  e,*  with  i*  =  argrnin^  Vh'(x)-  In 
short,  this  relaxation  is  not  exact  and  therefore  cannot  be  used  to  globally 
solve  the  original  problem. 
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5.2  Convex  relaxation  on  dual  variables  (Pock  et  al. 

in]) 

Another  decomposition  of  the  labeling  function  u:  Q  —>  (0, . . . ,  k}  is  the 
following.  Put  do  =  1,  dk+i  =  0,  and  for  i  —  1, . . . ,  k, 

«.w  =  {1  ifu(a:)-i’ 

I  0  otherwise. 

Then  u  =  JA  9%,  and  a  one-to-one  correspondence  between  the  labeling  func¬ 
tion  u  and  the  /c-tuple  (0l5 . . . ,  6k)  is  guaranteed  by  imposing  the  condition 
(9i,...,9k)  G  B0,  where 

B0  =  {(0i,...,0fc):  ^  -»•  {0,  l}fc :  0  <  9k(x)  <■■■<  9x{x)  <1  Va;  G  12}  . 

Note  the  relationship  between  this  decomposition  and  the  previous  one: 


9i(x)  -  9i+i(x)  =  if} i(x ). 


Under  this  framework,  the  segmentation  problem  is 

min  {V'  /  |V(0i(x)  -  0i+i(  ■*0)1  T  (9i(x)  9i-\-i(x)^)gi(x)  dx\  .  (19) 

(0i,...A)eso  J 

The  novel  approach  of  [29]  is  to  use  the  dual  formulation  of  TV  to  obtain 
a  more  tractable  minimization  problem  from  (19).  Recall  that  the  regular¬ 
ization  term  in  the  Mumford- Shah/ Potts  model  is  the  (d  —  l)-dimensional 
Hausdorff  measure  of  the  so-called  jump  set  Ju  (that  is,  the  set  of  essential 
discontinuities)  of  the  function  u  of  bounded  variation  (see  [17]  for  more  on 
BV  functions  and  their  properties).  In  other  words,  up  to  constant  factors, 


K  p 

£  /  |V(0i(x)-0i+1(a;))|  dx  =  Hd-\Ju). 

i= o  Jn 

Although  highly  nontrivial,  this  can  be  written  (see  the  technical  report  [11]) 

9i{x)  div^(x)  dx  \  , 


{k 

£-  [ 
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where 


ic  . &):<> 


' :  ^j(x)  <  1  \/x  G  fi,  1  <  ii  <  i2  <  k  >  . 

ii<i<j2  ) 


It  is  straightforward  to  see  that  this  set  /C  is  convex.  Thus,  this  yields  the 
problem 


min  sup  Wri  Oi(x^j  div T  {9i(pc)  dj+i^))^^)  , 

which  turns  out  to  be  equivalent  (again  see  [11])  to 


min  sup  <  -  /  9l{x)  div^(x)  +  \9i(x)  -  9i+1(x)\  g^x)  dx  >  , 

0eB  ^  t ~1  Jn 


where 


B  =  {{0u...,0k):n^{  0,  l}fc} 


is  binary  but  now  unordered.  The  problem  remains  non-convex  due  to  the 
set  B.  As  we  have  seen  previously,  a  relaxation  is  taken  to  minimize  over  9 
over  the  convex  hypercube 

R=  {(9u...,9k):  [0,l]fc}. 


However,  once  again  the  relaxation  fails  to  be  exact:  there  is  no  guarantee 
that  we  may  obtain  a  true  minimizer  over  B  from  a  minimizer  over  R. 


5.3  Exact,  separately  convex  relaxation 

Finally,  let  us  mention  one  more  alternative  approach  to  use  convex  relax¬ 
ation  for  the  multi-phase  problem.  Suppose  the  number  of  phases  k  +  1  is 
equal  to  2"  for  some  integer  n.  (In  fact,  this  is  no  restriction  whatsoever  since 
superfluous  phases  may  be  taken  to  be  empty.)  The  multi-phase  framework 
[34]  uses  n  level  set  functions  </q, . . . ,  tpn  to  represent  the  segmentation,  where 
the  union  of  the  zero  level-sets  of  the  </q  yields  the  edges  of  the  segmentation. 
Equivalently,  we  can  consider  the  binary  vector  of  Heaviside  functions  of  the 
(pi,  say  h  =  (hi, . . . ,  hn) :  0  — ^  {0,  l}n;  each  unique  binary  vector  corresponds 
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to  one  phase  of  the  partition.  The  segmentation  problem  may  be  formulated 
as 


(20) 


where 


^b(h)  =  Y[zbi(hi), 

i= 1 


ZbAv) 


1  -y  if  h  =  0, 
y  if  bi  =  1. 


The  region  descriptor  g b  is  simply  a  relabeling  of  gj  where  the  components 
of  b  are  the  base-2  digits  for  j .  For  example,  in  the  four-phase  case  when 
n  =  2, 


G(hi,h-2)  =  /  \Vhi\dx+  /  |  Vh2 1  dx  +  /  (x) hi (x) h2 (x)  dx 


'n 


+  /  g{Q^(x)(l^  h1(x))h2{x)dx+  /  g{lfi)(x)hi(x)(l  -  h2(x))  dx 
Jn  Jn 

+  /  9( o,o)(^)(l  -  /ii (a?))(l  -  h2(x))  dx. 

Jn 

ft  is  straightforward  to  prove  that  (20)  has  an  exact  convex  relaxation. 

Proposition  2.  We  have 

min  G(h)  =  min  G(h). 

he{o,i}n  he  [0,1]” 

Moreover,  if  h* :  hi  — >  [0,  l]n  is  a  solution  of  the  relaxed  problem,  then 

(l{/ll(x)>tl}>  ‘  ’  l{/l*  (x)>tn})  ‘  ^  >  {0) 

is  a  solution  of  the  binary  problem  (20)  for  any  (ti, . . . ,  tn)  G  (0,  l)n. 

Proof.  We  generalize  the  two-phase  proof  from  [7].  Given  a  solution  h*  of  the 
relaxed  problem,  define  <j)(x,t)  =  ti), . . . ,  (f>n(x,  tn)),  where  f>i(x,ti)  = 

1  {h*(x)>u}i  for  any  t  =  (G, . . .  ,tn)  e  [0,  l]n.  By  the  coarea  formula, 


\Vh*\  dx  = 


\V<f>i(x,t)\  dt. 


n  Jo 
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By  the  layer-cake  formula, 


Hence, 


Zb(h*)gb(x)dx 


'  [0,1]”  Jsi 


g\3(x)Z\3((j)(x,  t))  dx  dt. 


G(h*) 


G(4>(x,  t))  dt. 


Let  h  be  a  solution  of  the  binary  problem.  Since  (j)  €  {0,  l}n,  we  have 
G(4>)  >  G(h)  for  all  t.  On  the  other  hand,  because  the  relaxed  problem  is 
minimized  over  a  larger  set,  we  have  G(h)  >  G(h*).  Thus, 


G{ h*)=  f  G(<f>(x,t))dt  >  f  G{h)  =  G{h)>G{h*), 

which  implies  that  G(h)  =  G(h*)  =  G(4>(x,t))  for  all  t  e  (0,  l)n.  □ 

The  objective  functional  G(-)  is  not  convex,  but  it  is  separately  convex 
in  each  hi,  i.e. ,  for  all  i,  if  every  other  hj  ( j  ^  i)  is  fixed,  then  the  function 
is  convex  in  h% .  Separate  convexity  cannot  in  general  provide  an  algorithm 
that  guarantees  to  find  a  global  minimizer.  However,  in  practice,  this  method 
seems  to  give  satisfactory  results. 


5.4  Conclusion 

In  conclusion,  we  have  provided  a  method  that  guarantees  a  global  solution 
of  the  multi-phase  piecewise  constant  Mumford-Shah  segmentation  problem. 
The  method  solves  a  sequence  of  convex  optimization  problems  to  yield  a 
global  solution  of  the  original  non-convex  problem.  Our  method  appears 
to  be  the  first  in  the  literature  that  can  guarantee  that  a  global  solution  is 
obtained. 

In  this  section,  we  have  described  alternative  approaches  to  solving  this 
multi-phase  segmentation  problem.  In  some  sense,  compared  with  our  ap¬ 
proach,  these  methods  are  more  natural  multi-phase  generalizations  of  [12]. 
However,  each  of  these  methods  is  flawed  in  the  sense  that  a  global  minimizer 
cannot  be  guaranteed.  The  techniques  described  in  Sections  5. 1-5. 2  are  not 


23 


exact  relaxations;  indeed,  in  general,  one  cannot  obtain  a  global  minimizer 
of  the  multi-phase  problem  from  a  minimizer  of  the  respective  convex  refor¬ 
mulations.  On  the  other  hand,  the  technique  described  in  Section  5.3  is  an 
exact  relaxation,  but  unfortunately  the  lack  of  convexity  of  the  reformulation 
prohibits  the  assurance  of  a  global  minimizer. 

We  see  that  the  strength  of  our  approach  is  to  formulate  the  multi-phase 
problem  in  such  a  way  that  we  can  apply  a  convex  relaxation  approach  faith¬ 
fully.  This  involves  using  the  PCLSM  framework  with  the  augmented  La- 
grangian,  the  dual  formulation  of  TV,  and  functional  lifting  for  the  convex 
relaxation.  This  combination  of  tools  is  able  to  guarantee  the  global  solution, 
in  contrast  to  previous  methods. 

There  is  still  future  work  that  can  be  done  in  this  direction.  This  includes 
developing  algorithms  to  improve  the  efficiency  of  our  method.  We  also  want 
to  consider  the  more  difficult  problem  in  which  the  optimal  constants  (i.e., 
the  region  descriptors)  are  not  known  a  priori. 
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